Figure 1 - This figures shows…
library(png)
library(maptools)
Checking rgeos availability: TRUE
library(raster)
library(gam)
Loading required package: splines
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded gam 1.14
setwd("~/Desktop/Botero postdoc 2016/Human density and the origins of agriculture/")
par(mar=c(0,0,0,20))
d <- readPNG("Larson_dates.png")
plot(seq(0,18, length.out = 19), seq(0,36, length.out = 19), type="n",ylim=c(0,36),xlim=c(0, 18), xaxt="n")
rasterImage(d, 0,0,18,36, interpolate=TRUE, col=d)
Start_of_early_window <- 16-12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 17-4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
These dates are provided in the supplimentary information for the Larson (2014) paper. I’ve copied those values into a .csv table provided here.
domestication_times <- read.csv("Domestication timing larson 2014.csv")
dim(domestication_times)
[1] 77 8
| Region | Species | Start.Exploitation | Finish.Exploitation | Start.predomestication | Finish.predomestication | Start.Domestication | Finish.Domestication |
|---|---|---|---|---|---|---|---|
| Southwest asia | Wheat | 12.00 | 11.25 | 11.25 | 11.00 | 11.00 | 9.00 |
| Southwest asia | Barley | 12.00 | 11.25 | 11.25 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Lentil | 12.00 | 11.00 | 11.00 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Pea | 11.50 | 11.00 | 11.00 | 10.00 | 10.00 | 8.50 |
| Southwest asia | Chickpea | 11.00 | 10.50 | 10.50 | 10.25 | 10.25 | 8.25 |
| Southwest asia | Broadbean | NA | NA | NA | NA | 10.50 | NA |
| Southwest asia | Flax | 12.00 | 9.50 | NA | NA | 9.50 | NA |
| Southwest asia | Olive | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| Southwest asia | Sheep | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Goat | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Pig | 12.00 | 11.50 | 11.50 | 9.75 | 10.25 | 9.00 |
| Southwest asia | Cattle, taurine | 11.50 | 10.50 | 10.50 | 10.25 | 10.25 | 8.00 |
| Southwest asia | Cat | NA | NA | 10.50 | 4.00 | 4.00 | NA |
| South Asia | Tree cotton | 8.50 | 4.50 | NA | NA | 4.50 | NA |
| South Asia | Rice | 8.00 | 5.00 | 5.00 | 4.00 | 4.00 | 2.50 |
| South Asia | Little millet | NA | NA | NA | NA | 4.50 | NA |
| South Asia | Browntop millet | NA | NA | NA | NA | 4.00 | NA |
| South Asia | Mungbean | NA | NA | 4.50 | 3.50 | 3.50 | 3.00 |
| South Asia | Pigeonpea | NA | NA | NA | NA | 3.50 | NA |
| South Asia | Zebu cattle | 9.00 | 8.00 | NA | NA | 8.00 | 6.50 |
| South Asia | Water buffalo | 6.00 | 4.50 | NA | NA | 4.50 | NA |
| East Asia | Broomcorn Millet | 10.00 | 8.00 | NA | NA | 8.00 | NA |
| East Asia | Foxtail millet | 11.50 | 7.50 | NA | NA | 7.50 | NA |
| East Asia | Rice | 10.00 | 8.00 | 8.00 | 7.50 | 7.50 | 5.00 |
| East Asia | Soybean | 8.50 | 5.50 | NA | NA | 5.50 | 4.00 |
| East Asia | Ramie | NA | NA | NA | NA | 5.25 | NA |
| East Asia | Melon | 7.00 | 4.00 | NA | NA | 4.00 | 3.75 |
| East Asia | Pig | 12.00 | 8.50 | NA | NA | 8.50 | 6.00 |
| East Asia | Silkworm | 7.00 | 5.25 | NA | NA | 5.25 | NA |
| East Asia | Yak | NA | NA | NA | NA | 4.25 | NA |
| East Asia | Horse | 7.50 | 6.75 | 6.75 | 5.50 | 5.50 | 4.00 |
| East Asia | Bactrian Camel | NA | NA | NA | NA | 4.50 | NA |
| East Asia | Duck | 2.50 | 1.00 | NA | NA | 1.00 | NA |
| East Asia | Chicken | 6.00 | 4.00 | NA | NA | 4.00 | NA |
| New Guinea | Banana | 10.00 | 7.00 | 7.00 | 4.00 | 4.00 | NA |
| New Guinea | Taro | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| New Guinea | Yam | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| Africa and Arabia | Date palm | 7.00 | 6.00 | NA | NA | 5.00 | NA |
| Africa and Arabia | Sorghum | 8.00 | 4.00 | NA | NA | 4.00 | NA |
| Africa and Arabia | Pearl millet | NA | NA | NA | NA | 4.50 | 3.50 |
| Africa and Arabia | Fonio | NA | NA | NA | NA | 2.50 | NA |
| Africa and Arabia | Cowpea | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Hyacinth bean | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Rice | 3.50 | 2.00 | NA | NA | 2.00 | NA |
| Africa and Arabia | Oil palm | 9.25 | 3.50 | NA | NA | 3.50 | NA |
| Africa and Arabia | Cattle | NA | NA | 9.00 | 7.75 | 7.75 | 6.50 |
| Africa and Arabia | Donkey | 9.00 | 5.50 | NA | NA | 5.50 | 3.50 |
| Africa and Arabia | Dromedary camel | 6.50 | 3.00 | NA | NA | 3.00 | NA |
| Africa and Arabia | Guinea fowl | NA | NA | 2.50 | 1.50 | 1.50 | NA |
| North America | Squash | 6.50 | 5.00 | NA | NA | 5.00 | NA |
| North America | Sunflower | 6.00 | 4.75 | NA | NA | 4.00 | NA |
| North America | Sumpweed | 6.00 | 4.50 | NA | NA | 4.00 | NA |
| North America | Pitseed goosefoot | 4.75 | 3.75 | NA | NA | 3.75 | NA |
| Meso-america | Squash (pepo) | NA | NA | NA | NA | 10.00 | 9.50 |
| Meso-america | Maize | 10.00 | 9.00 | NA | NA | 9.00 | NA |
| Meso-america | Foxtail millet-grass | NA | NA | NA | NA | 6.00 | 4.00 |
| Meso-america | Common bean | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Avocado | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Chile pepper | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Turkey | NA | NA | NA | NA | 2.00 | NA |
| South America | Chili pepper | NA | NA | NA | NA | 6.00 | NA |
| South America | Peanut | NA | NA | NA | NA | 5.00 | NA |
| South America | Cotton | NA | NA | NA | NA | 6.00 | NA |
| South America | Coca | NA | NA | NA | NA | 8.00 | NA |
| South America | Now-minor root crops (arrowroot, leren) | NA | NA | NA | NA | 9.00 | NA |
| South America | Squash | NA | NA | NA | NA | 10.00 | NA |
| South America | Common bean | NA | NA | NA | NA | 5.00 | NA |
| South America | Lima bean | NA | NA | 8.25 | NA | 6.00 | NA |
| South America | Monioc | NA | NA | NA | NA | 7.00 | NA |
| South America | Sweet potato | NA | NA | NA | NA | 5.00 | NA |
| South America | White potato | 7.00 | 4.50 | NA | NA | 4.00 | NA |
| South America | Quinoa | 5.00 | NA | NA | NA | 3.50 | NA |
| South America | Yam | NA | NA | NA | NA | 5.50 | NA |
| South America | Llama | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| South America | Alpaca | 10.00 | 5.00 | NA | NA | 5.00 | NA |
| South America | Guinea pig | NA | NA | NA | NA | 5.00 | NA |
| South America | Muscovy Duck | NA | NA | NA | NA | 4.00 | NA |
par(mar=c(5,4,6,1))
dates <- unlist(domestication_times[3:8])
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset" )
mtext("This tells us about how evenly our evidence is distributed in time", 3, line=1)
hist(dates, breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset with Larson(2014) date windows")
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
mtext("Early Holocene", 3, line = -1, adj=.3)
mtext("Middle Holocene", 3, line= -1, adj=.6)
par(mfrow=c(2,3), mar=c(4,4,2,0))
dim(domestication_times)
[1] 77 8
specific_dates <- domestication_times[,3:8]
for(i in c(1, 3, 5, 2, 4, 6)){
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main= names(specific_dates)[i])
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
}
I’m creating new rows for this table, combining dates in different ways to make the CDFs below look more authentic. This makes it so that pre-ag always happens before post-ag. What I’ve done is given the later date to the earlier date when those dates are missing.
h <- which(is.na(domestication_times[,3]))
domestication_times <- cbind(domestication_times, rep(NA, length(domestication_times[,1])))
domestication_times[,9] <- domestication_times[,3]
domestication_times[h,9] <- domestication_times[h,7]
colnames(domestication_times)[9] <- "adopt exploitation date"
domestication_times[,10] <- domestication_times[,7]
domestication_times[which(is.na(domestication_times[,10])),10] <- 0
colnames(domestication_times)[10] <- "start of ag"
#save(domestication_times, file="~/Desktop/Human density and the origins of agriculture/Domestication timing larson 2014.Rdata")
I think these are best described by a cummulative distribution, showing how they accumulate over time.
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(levels(domestication_times$Region)[ type_number])
print(match)
print(j)
}
[1] "Africa and Arabia"
[1] 7.00 8.00 4.50 2.50 3.75 3.75 3.50 9.25 7.75 9.00 6.50 1.50
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
[1] "East Asia"
[1] 10.00 11.50 10.00 8.50 5.25 7.00 12.00 7.00 4.25 7.50 4.50 2.50 6.00
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
[1] "Meso-america"
[1] 10 10 6 3 3 3 2
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
[1] "New Guinea"
[1] 10 10 10
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
[1] "North America"
[1] 6.50 6.00 6.00 4.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
[1] "South America"
[1] 6.0 5.0 6.0 8.0 9.0 10.0 5.0 6.0 7.0 5.0 7.0 5.0 5.5 10.0 10.0 5.0 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
[1] "South Asia"
[1] 8.5 8.0 4.5 4.0 3.5 3.5 9.0 6.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
[1] "Southwest asia"
[1] 12.0 12.0 12.0 11.5 11.0 10.5 12.0 10.0 12.0 12.0 12.0 11.5 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
if(i == 2 | i == 1)axis(2)
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
}
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
abline(v= maxer - quantile(j)[2], col="limegreen", lwd=2)
if(i == 2 | i == 1)axis(2)
if(i == 2)mtext("25%", 3, line=3.5, adj=-1, col="limegreen")
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
if(i == 4)mtext("Choose a y to predict an x", 3, line=3.3, col="cornflowerblue")
break_one <- maxer
break_two <- maxer - quantile(j)[2]
polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.2), border=adjustcolor("cornflowerblue",alpha=1))
lines(x=c(break_two, break_two), y=c(0,-1), col="cornflowerblue")
abline(h = 25, col="limegreen", lwd=2)
}
Make this a function. There is a choice of two methods here. At the end of this section we need to print the desision we’re passing to the later analyses.
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)
[1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" NA "E_North_Ameri"
[7] "New_Guinea" "Mesoamerica" "N_Lowland_SA" "NW_Lowland_SA" "Sava_W_India" "S_India"
[13] "Ganges_E_Indi" "Chinese_loess" "Japanese" "Lower-MiddleY" "South trop ch" NA
[19] "Southwes amaz" "C/S_Andes"
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA <NA> <NA> New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
18 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Ganges_E_Indi ... West Africa T
origins_subset$name
NULL
library(maps)
map()
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
map()
d <- readPNG("Larson_origins.png")
rasterImage(d, -180, -90, 180, 110, interpolate=TRUE, col=d)
map(add=TRUE)
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
# need to reproject
This is obviously a bad projection fit right now.
#subset and reorder origins. This is currently done at the end of the plot but should be moved forward.
# Load data for population density
load("PopD_all_December.rdata")
PopD.ALL
class : RasterStack
dimensions : 288, 720, 207360, 18 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : fourK, fiveK, sixK, sevenK, eightK, nineK, tenK, elevenK, twelveK, thirteenK, fourteenK, fifteenK, sixteenK, seventeenK, eighteenK, ...
min values : 5.611358e-07, 1.067142e-06, 2.508241e-06, 6.317553e-06, 2.286934e-05, 7.631922e-05, 1.272693e-04, 2.118215e-04, 2.602175e-04, 3.226203e-04, 4.390267e-04, 5.572032e-04, 7.313966e-04, 8.286005e-04, 8.297062e-04, ...
max values : 2.051069, 2.013452, 2.142908, 1.888403, 1.863014, 1.880628, 1.650615, 1.678033, 1.697732, 1.499115, 1.517264, 1.443677, 1.464867, 1.453581, 1.436394, ...
# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
r
class : RasterLayer
dimensions : 288, 720, 207360 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : fourK
values : 5.611358e-07, 2.051069 (min, max)
We need to justify our decision to use a GAM over other models. This should include citations to back up those arguments.
We should make our decisions very transparent here. We should be able to justify our decision of 3 degrees of freedom over other possible values.
# Read the polygons
origins <- readShapePoly('Origins_updated.shp')
# Extract data
per.origin <- extract(r, origins, cellnumber = TRUE, buffer = 100000)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
names(per.origin) <- origins@data[, 1]
str(per.origin)
List of 20
$ W_African_Sav: num [1:309, 1:2] 92506 92507 92508 92509 92510 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Sudanic_Savan: num [1:306, 1:2] 99050 99051 99052 99053 99054 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ West Africa T: num [1:427, 1:2] 102609 102610 102611 103326 103327 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Ethipian plat: num [1:275, 1:2] 106281 106282 106283 106998 106999 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NA : num [1:195, 1:2] 67404 67405 67406 68119 68120 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ E_North_Ameri: num [1:170, 1:2] 64270 64271 64984 64985 64986 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ New_Guinea : num [1:24, 1:2] 127360 127361 127362 128080 128081 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Mesoamerica : num [1:73, 1:2] 93032 93033 93034 93035 93036 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ N_Lowland_SA : num [1:39, 1:2] 110373 110374 110375 111092 111093 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NW_Lowland_SA: num [1:24, 1:2] 123319 123320 123321 123322 123323 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Sava_W_India : num [1:34, 1:2] 83312 84031 84032 84749 84750 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ S_India : num [1:18, 1:2] 99153 99154 99872 99873 99874 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Ganges_E_Indi: num [1:92, 1:2] 85494 85495 85496 85497 85498 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Chinese_loess: num [1:84, 1:2] 72598 72599 73318 73319 73320 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Japanese : num [1:36, 1:2] 59681 59682 60401 61121 61843 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Lower-MiddleY: num [1:131, 1:2] 72547 72548 72549 73267 73268 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ South trop ch: num [1:178, 1:2] 84818 84819 84820 84821 84822 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NA : num [1:258, 1:2] 62493 62494 62495 62496 62497 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Southwes amaz: num [1:194, 1:2] 137758 137759 137760 137761 137762 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ C/S_Andes : num [1:165, 1:2] 137736 137737 137738 137739 137740 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
all_vectors <- rep(NA, 3)
for(h in 1:20){
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- as.vector(originI)
all_vectors_pre <- cbind(rep(names(per.origin)[h], length(x_value_vector)),x_value_vector, y_value_vector)
all_vectors <- rbind(all_vectors, all_vectors_pre)
}
all_vectors <- as.data.frame(all_vectors)
all_vectors <- cbind(all_vectors, scale(as.numeric(all_vectors[,3])))
some row.names duplicated: 3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,336,337,338,339,340,341,342,343,344,345,346,347,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,367,368,369,370,371,372,373,374,375,376,377,378,379,380,381,382,383,384,385,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420,421,422,423,424,425,426,427,428,429,430,431,432,433,434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,450,451,452,453,454,455,456,457,458,459,460,461,462,463,464,465,466,467,468,469,470,471,472,473,474,475,476,477,478,479,480,481,482,483,484,485,486,487,488,489,490,491,492,493,494,495,496,497,498,499,500,501,502,503,504,505,506,507,508,509,510,511,512,513,514,515,516,517,518,519,520,521,522,523,524,525,526,527,528,529,530,531,532,533,534,535,536,537,538,539,540,541,542,543,544,545,546,547,548,549,550,551,552,553,554,555,556,557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606,607,608,609,610,611,612,613,614,615,616,617,618,619,620,621,622,623,624,625,626,627,628,629,630,631,632,633,634,635,636,637,638,639,640,641,642,643,644,645,646,647,648,649,650,651,652,653,654,655,656,657,658,659,660,661,662,663,664,665,666,667,668,669,670,671,672,673,674,675,676,677,678,679,680,681,682,683,684,685,686,687,688,689,690,691,692,693,694,695,696,697,698,699,700,701,702,703,704,705,706,707,708,709,710,711,712,713,714,715,716,717,718,719,720,721,722,723,724,725,726,727,728,729,730,731,732,733,734,735,736,737,738,739,740,741,742,743,744,745,746,747,748,749,750,751,752,753,754,755,756,757,758,759,760,761,762,763,764,765,766,767,768,769,770,771,772,773,774,775,776,777,778,779,780,781,782,783,784,785,786,787,788,789,790,791,792,793,794,795,796,797,798,799,800,801,802,803,804,805,806,807,808,809,810,811,812,813,814,815,816,817,818,819,820,821,822,823,824,825,826,827,828,829,830,831,832,833,834,835,836,837,838,839,840,841,842,843,844,845,846,847,848,849,850,851,852,853,854,855,856,857,858,859,860,861,862,863,864,865,866,867,868,869,870,871,872,873,874,875,876,877,878,879,880,881,882,883,884,885,886,887,888,889,890,891,892,893,894,895,896,897,898,899,900,901,902,903,904,905,906,907,908,909,910,911,912,913,914,915,916,917,918,919,920,921,922,923,924,925,926,927,928,929,930,931,932,933,934,935,936,937,938,939,940,941,942,943,944,945,946,947,948,949,950,951,952,953,954,955,956,957,958,959,960,961,962,963,964,965,966,967,968,969,970,971,972,973,974,975,976,977,978,979,980,981,982,983,984,985,986,987,988,989,990,991,992,993,994,995,996,997,998,999,1000,1001,1002,1003,1004,1005,1006,1007,1008,1009,1010,1011,1012,1013,1014,1015,1016,1017,1018,1019,1020,1021,1022,1023,1024,1025,1026,1027,1028,1029,1030,1031,1032,1033,1034,1035,1036,1037,1038,1039,1040,1041,1042,1043,1044,1045,1046,1047,1048,1049,1050,1051,1052,1053,1054,1055,1056,1057,1058,1059,1060,1061,1062,1063,1064,1065,1066,1067,1068,1069,1070,1071,1072,1073,1074,1075,1076,1077,1078,1079,1080,1081,1082,1083,1084,1085,1086,1087,1088,1089,1090,1091,1092,1093,1094,1095,1096,1097,1098,1099,1100,1101,1102,1103,1104,1105,1106,1107,1108,1109,1110,1111,1112,1113,1114,1115,1116,1117,1118,1119,1120,1121,1122,1123,1124,1125,1126,1127,1128,1129,1130,1131,1132,1133,1134,1135,1136,1137,1138,1139,1140,1141,1142,1143,1144,1145,1146,1147,1148,1149,1150,1151,1152,1153,1154,1155,1156,1157,1158,1159,1160,1161,1162,1163,1164,1165,1166,1167,1168,1169,1170,1171,1172,1173,1174,1175,1176,1177,1178,1179,1180,1181,1182,1183,1184,1185,1186,1187,1188,1189,1190,1191,1192,1193,1194,1195,1196,1197,1198,1199,1200,1201,1202,1203,1204,1205,1206,1207,1208,1209,1210,1211,1212,1213,1214,1215,1216,1217,1218,1219,1220,1221,1222,1223,1224,1225,1226,1227,1228,1229,1230,1231,1232,1233,1234,1235,1236,1237,1238,1239,1240,1241,1242,1243,1244,1245,1246,1247,1248,1249,1250,1251,1252,1253,1254,1255,1256,1257,1258,1259,1260,1261,1262,1263,1264,1265,1266,1267,1268,1269,1270,1271,1272,1273,1274,1275,1276,1277,1278,1279,1280,1281,1282,1283,1284,1285,1286,1287,1288,1289,1290,1291,1292,1293,1294,1295,1296,1297,1298,1299,1300,1301,1302,1303,1304,1305,1306,1307,1308,1309,1310,1311,1312,1313,1314,1315,1316,1317,1318,1319,1320,1321,1322,1323,1324,1325,1326,1327,1328,1329,1330,1331,1332,1333,1334,1335,1336,1337,1338,1339,1340,1341,1342,1343,1344,1345,1346,1347,1348,1349,1350,1351,1352,1353,1354,1355,1356,1357,1358,1359,1360,1361,1362,1363,1364,1365,1366,1367,1368,1369,1370,1371,1372,1373,1374,1375,1376,1377,1378,1379,1380,1381,1382,1383,1384,1385,1386,1387,1388,1389,1390,1391,1392,1393,1394,1395,1396,1397,1398,1399,1400,1401,1402,1403,1404,1405,1406,1407,1408,1409,1410,1411,1412,1413,1414,1415,1416,1417,1418,1419,1420,1421,1422,1423,1424,1425,1426,1427,1428,1429,1430,1431,1432,1433,1434,1435,1436,1437,1438,1439,1440,1441,1442,1443,1444,1445,1446,1447,1448,1449,1450,1451,1452,1453,1454,1455,1456,1457,1458,1459,1460,1461,1462,1463,1464,1465,1466,1467,1468,1469,1470,1471,1472,1473,1474,1475,1476,1477,1478,1479,1480,1481,1482,1483,1484,1485,1486,1487,1488,1489,1490,1491,1492,1493,1494,1495,1496,1497,1498,1499,1500,1501,1502,1503,1504,1505,1506,1507,1508,1509,1510,1511,1512,1513,1514,1515,1516,1517,1518,1519,1520,1521,1522,1523,1524,1525,1526,1527,1528,1529,1530,1531,1532,1533,1534,1535,1536,1537,1538,1539,1540,1541,1542,1543,1544,1545,1546,1547,1548,1549,1550,1551,1552,1553,1554,1555,1556,1557,1558,1559,1560,1561,1562,1563,1564,1565,1566,1567,1568,1569,1570,1571,1572,1573,1574,1575,1576,1577,1578,1579,1580,1581,1582,1583,1584,1585,1586,1587,1588,1589,1590,1591,1592,1593,1594,1595,1596,1597,1598,1599,1600,1601,1602,1603,1604,1605,1606,1607,1608,1609,1610,1611,1612,1613,1614,1615,1616,1617,1618,1619,1620,1621,1622,1623,1624,1625,1626,1627,1628,1629,1630,1631,1632,1633,1634,1635,1636,1637,1638,1639,1640,1641,1642,1643,1644,1645,1646,1647,1648,1649,1650,1651,1652,1653,1654,1655,1656,1657,1658,1659,1660,1661,1662,1663,1664,1665,1666,1667,1668,1669,1670,1671,1672,1673,1674,1675,1676,1677,1678,1679,1680,1681,1682,1683,1684,1685,1686,1687,1688,1689,1690,1691,1692,1693,1694,1695,1696,1697,1698,1699,1700,1701,1702,1703,1704,1705,1706,1707,1708,1709,1710,1711,1712,1713,1714,1715,1716,1717,1718,1719,1720,1721,1722,1723,1724,1725,1726,1727,1728,1729,1730,1731,1732,1733,1734,1735,1736,1737,1738,1739,1740,1741,1742,1743,1744,1745,1746,1747,1748,1749,1750,1751,1752,1753,1754,1755,1756,1757,1758,1759,1760,1761,1762,1763,1764,1765,1766,1767,1768,1769,1770,1771,1772,1773,1774,1775,1776,1777,1778,1779,1780,1781,1782,1783,1784,1785,1786,1787,1788,1789,1790,1791,1792,1793,1794,1795,1796,1797,1798,1799,1800,1801,1802,1803,1804,1805,1806,1807,1808,1809,1810,1811,1812,1813,1814,1815,1816,1817,1818,1819,1820,1821,1822,1823,1824,1825,1826,1827,1828,1829,1830,1831,1832,1833,1834,1835,1836,1837,1838,1839,1840,1841,1842,1843,1844,1845,1846,1847,1848,1849,1850,1851,1852,1853,1854,1855
colnames(all_vectors) <- c("location_name", "x_values", "density_values", "scaled_density_values")
all_vectors
not.origin <- r[r != per.origin[[1]][,1]]
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
longer object length is not a multiple of shorter object length
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:18){
#h <- 3
density_trend <- all_vectors[which(all_vectors[,1] == levels(all_vectors[,1])[h]),]
plot(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", ylim=c(-5,5), xlim=c(4,21), type="n")
points(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), cex=0.5, col=adjustcolor("cornflowerblue", 0.5))
}
#all origins
plot(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", ylim=c(-5,5), xlim=c(4,21), type="n")
points(as.numeric(all_vectors[,2]), as.numeric(all_vectors[,4]), cex=0.5, col=adjustcolor("cornflowerblue", 0.5))
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:18){
#h <- 3
density_trend <- all_vectors[which(all_vectors[,1] == levels(all_vectors[,1])[h]),]
plot(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", ylim=c(-3,3), xlim=c(4,21), type="n", xaxt="n", yaxt="n")
points(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), cex=0.5, col=adjustcolor("grey", 0.5))
ordered_density_trend <- density_trend[order(density_trend$x_values), ]
gammer <- loess(as.numeric(ordered_density_trend[,4]) ~ as.numeric(ordered_density_trend[,2]))
summary(gammer)
lines(as.numeric(ordered_density_trend[,2]) ,predict(gammer), col="cornflowerblue", lwd=2)
}
density_trend <- all_vectors[-1,]
plot(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", ylim=c(-3,3), xlim=c(4,21), type="n", xaxt="n", yaxt="n")
points(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), cex=0.5, col=adjustcolor("grey", 0.5))
ordered_density_trend <- density_trend[order(density_trend$x_values), ]
gammer <- loess(as.numeric(ordered_density_trend[,4]) ~ as.numeric(ordered_density_trend[,2]))
summary(gammer)
Call:
loess(formula = as.numeric(ordered_density_trend[, 4]) ~ as.numeric(ordered_density_trend[,
2]))
Number of Observations: 54576
Equivalent Number of Parameters: 4.29
Residual Standard Error: 0.9855
Trace of smoother matrix: 4.66 (exact)
Control settings:
span : 0.75
degree : 2
family : gaussian
surface : interpolate cell = 0.2
normalize: TRUE
parametric: FALSE
drop.square: FALSE
length(ordered_density_trend[,2])
[1] 54576
length(predict(gammer))
[1] 54576
lines(as.numeric(ordered_density_trend[,2]) ,predict(gammer), col="cornflowerblue", lwd=2)
span
the parameter α which controls the degree of smoothing.
The size of the neighbourhood is controlled by α (set by span or enp.target). For α < 1, the neighbourhood includes proportion α of the points, and these have tricubic weighting (proportional to (1 - (dist/maxdist)3)3). For α > 1, all points are used, with the ‘maximum distance’ assumed to be α^(1/p) times the actual maximum distance for p explanatory variables.
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:18){
#h <- 3
density_trend <- all_vectors[which(all_vectors[,1] == levels(all_vectors[,1])[h]),]
ordered_density_trend <- density_trend[order(density_trend$x_values), ]
gammer <- loess(as.numeric(ordered_density_trend[,4]) ~ as.numeric(ordered_density_trend[,2]), span = .5)
summary(gammer)
predict_gam <- predict(gammer, se=TRUE)
plot(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", ylim=c(-3,3), xlim=c(0,21), type="n", xaxt="n", yaxt="n")
#axis(1, label= seq(4,21, by=1), at=rev(seq(1,18, by=1)))
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
abline(h=0, lty=2, col="grey")
lines(ordered_density_trend[,2] ,predict_gam$fit, col="cornflowerblue")
lines(ordered_density_trend[,2] , predict_gam$fit + predict_gam$se.fit, col="grey", lty=1)
lines(ordered_density_trend[,2] , predict_gam$fit - predict_gam$se.fit, col="grey", lty=1)
mtext(as.character(ordered_density_trend[1,1]), 3, line=-1, cex=0.5)
}
par(mfrow=c(2,5), mar=c(0,0,0,0))
for(j in seq(0.2, 2, by=.2)){
plot(as.numeric(density_trend[,2]), as.numeric(density_trend[,4]), col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", ylim=c(-3,3), xlim=c(0,21), type="n", xaxt="n")
axis(1, label= seq(4,21, by=1), at=rev(seq(1,18, by=1)))
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
for(h in 1:18){
#h <- 3
density_trend <- all_vectors[which(all_vectors[,1] == levels(all_vectors[,1])[h]),]
ordered_density_trend <- density_trend[order(density_trend$x_values), ]
gammer <- loess(as.numeric(ordered_density_trend[,4]) ~ as.numeric(ordered_density_trend[,2]), span = j)
summary(gammer)
predict_gam <- predict(gammer, se=TRUE)
abline(h=0, lty=2, col="grey")
lines(ordered_density_trend[,2] ,predict_gam$fit, col="cornflowerblue")
}
mtext(paste("span = ", j), 3, line=-3, cex=0.5)
}
# need to add a global mean, an everything but the origins mean, and a buffer around the origins mean.
# Function standardization
std <- function(x) {
b <- (x - min(x)) / (max(x) - min(x))
return(rev(b))
}
diff_df <- function(h){
# Calculating mean and
global.means <- global.SD <- list()
for (j in 1:length(per.origin)) {
#print(j)
originI <- Pop[per.origin[[j]][, 1], ]
time <- 21:4
originI <- na.exclude(originI)
b <- apply(originI, 1, std)
nJ <- nrow(originI)
predictions <- matrix(nrow = nJ, ncol = length(time))
colnames(predictions) <- as.character(time)
for(i in 1:nJ) {
# Need to show a gradient of these df values.
model <- gam(b[, i] ~ s(time, df = h))
col <- sample(rainbow(100), 1)
predictions[i, ] <- predict(model)
}
global.means[[j]] <- apply(predictions, 2, mean)
global.SD[[j]] <- apply(predictions, 2, sd)
}
names(global.means) <- paste(names(per.origin), "Means")
names(global.SD) <- paste(names(per.origin), "SD")
return(list(global.means, global.SD))
}
# Confirm that the x-axis is oriented correctly
for_3 <- diff_df(1)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
i <- 1
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,20))
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
originI <- Pop[per.origin[[1]][, 1], ]
#colnames(originI)
#plot(21:4, rep(0, length(21:4)), type="n", ylim=c(0,1), xlim=c(21,4), xaxt="n", yaxt="n")
for(j in 1:18){
points(23- c(jitter(rep(4, length(originI[,j])), 4.5) + j), originI[,j], pch=19, cex=.1)
}
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]], type="b", pch=names(global.means[[i]]))
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[1], 3, line=-2)
means_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(means_matrix) <- c("origin", rev(seq(4, 21, by=1)))
means_matrix[,1] <- names(global.means)
for(i in 1:20){
means_matrix[i,2:19] <- global.means[[i]]
}
kable(means_matrix, caption= "Mean values")
| origin | 21 | 20 | 19 | 18 | 17 | 16 | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| W_African_Sav Means | 0.275750492966357 | 0.237971703906143 | 0.200504339081811 | 0.165047402018049 | 0.135516039740239 | 0.121748343406075 | 0.132429431895417 | 0.17103635569842 | 0.235948736377181 | 0.321617048802923 | 0.418964971526823 | 0.514676226549989 | 0.591121916502731 | 0.628378724452797 | 0.621163452983342 | 0.582047087468811 | 0.532681459452571 | 0.482909221409947 |
| Sudanic_Savan Means | -0.0186911478417769 | -0.00828615401443568 | 0.00413446846041128 | 0.0215161762523414 | 0.0473073567212149 | 0.0847812201109172 | 0.136209309975754 | 0.20228390605606 | 0.282347797217804 | 0.374044632344897 | 0.472296827728429 | 0.568020597195695 | 0.648225465565203 | 0.697539021870325 | 0.713328901689326 | 0.707327373243822 | 0.697887725176853 | 0.69056153763021 |
| West Africa T Means | -0.0145401490853128 | -0.00333603963380315 | 0.00977013345024412 | 0.027586956435481 | 0.0533165943413057 | 0.0895594750479319 | 0.138174778488521 | 0.199879802842878 | 0.274517388507917 | 0.360833593695506 | 0.45568720488448 | 0.552897898880818 | 0.643252356019774 | 0.715741781485123 | 0.768806992435192 | 0.809952043255802 | 0.850519193881187 | 0.893430170662931 |
| Ethipian plat Means | 0.0516408507031156 | 0.0622592531360626 | 0.0750787355293174 | 0.0933555709623621 | 0.12092332890413 | 0.16140265228906 | 0.217329358000213 | 0.28895304336974 | 0.373944304902601 | 0.467450912361356 | 0.562582640478741 | 0.650493007468181 | 0.722074582042064 | 0.768618080463592 | 0.790630347643354 | 0.797506826970228 | 0.803240018411344 | 0.812789660231236 |
| NA Means | 0.72808776413489 | 0.6723269776978 | 0.613165929035697 | 0.547108730705425 | 0.472515896153583 | 0.393590592158964 | 0.319669766181732 | 0.262708374009385 | 0.233548745835161 | 0.239616095881227 | 0.28323109870948 | 0.359870592475047 | 0.454992276384195 | 0.545528134731924 | 0.616545768777106 | 0.664148269449493 | 0.691583655592282 | 0.708950521834131 |
| E_North_Ameri Means | -0.0424144312128647 | -0.0202618739720356 | 0.00475194649595462 | 0.0365915191653316 | 0.079541401163135 | 0.136467523436837 | 0.207873937132646 | 0.291025513809716 | 0.380397583621781 | 0.468564149521066 | 0.547361392244133 | 0.609663875755146 | 0.654835785797328 | 0.689283124023394 | 0.721880145858871 | 0.761679332218606 | 0.814909663949399 | 0.875153810181115 |
| New_Guinea Means | 0.304105470512681 | 0.28132745701113 | 0.25734789630324 | 0.230996250201583 | 0.201793192056862 | 0.171544911503563 | 0.144548456661618 | 0.127184052406642 | 0.126877035836262 | 0.15052865629083 | 0.203062414319868 | 0.285273581529037 | 0.389014353820968 | 0.496874952215665 | 0.597210277530265 | 0.683771044637899 | 0.754313180721177 | 0.816774228406938 |
| Mesoamerica Means | -0.0315889824677606 | -0.00787599693194719 | 0.0184956204628149 | 0.0511692835844435 | 0.0938250236223921 | 0.148344738459587 | 0.214601573425648 | 0.290380258066549 | 0.372334867390379 | 0.456479966860125 | 0.538354702268212 | 0.613001244481024 | 0.676612820283987 | 0.726625668176205 | 0.767020170677251 | 0.80679724213028 | 0.855043623685109 | 0.908351302981932 |
| N_Lowland_SA Means | 0.208792535197897 | 0.197230324856633 | 0.185079121786267 | 0.171696272010441 | 0.156804931838545 | 0.141342946965234 | 0.128176787393512 | 0.122340507893562 | 0.130435008736839 | 0.159548596595071 | 0.215485053059509 | 0.300036154887671 | 0.405861213383308 | 0.51607541749585 | 0.62000422933541 | 0.712583747106496 | 0.792272847070064 | 0.866211062729842 |
| NW_Lowland_SA Means | 0.326576386801442 | 0.312436135756889 | 0.296990439986673 | 0.278723633997093 | 0.256416348833303 | 0.230169366346516 | 0.201945873000531 | 0.175714919736147 | 0.157651360929879 | 0.155352073191265 | 0.176275702224548 | 0.225642677767534 | 0.302153168560349 | 0.396329228756827 | 0.497310353388214 | 0.594432476311419 | 0.67925446713634 | 0.756553171147025 |
| Sava_W_India Means | -0.0237936944828927 | -0.00818464775151594 | 0.00933560457136895 | 0.0313826871867378 | 0.0607897006722834 | 0.0992817835941126 | 0.147568966381951 | 0.205468601119744 | 0.272443196878899 | 0.347413657258897 | 0.428584371379107 | 0.513039004997093 | 0.595722652725526 | 0.669265301648986 | 0.731180489316986 | 0.788223314024364 | 0.850666362502511 | 0.917274849593167 |
| S_India Means | -0.0222028081854339 | -0.0066002592031116 | 0.0109240264309662 | 0.0332350046742644 | 0.0636053304857559 | 0.104248025366543 | 0.156161312176903 | 0.218860716380081 | 0.291071667647164 | 0.371586637639489 | 0.458699358753368 | 0.549398894270441 | 0.638133144838091 | 0.71736832702396 | 0.786036841922908 | 0.849257343160103 | 0.913846434565112 | 0.979674443478567 |
| Ganges_E_Indi Means | -0.0168749067644347 | -0.00665207444708958 | 0.00516238304586039 | 0.0208799296723999 | 0.0431774421109922 | 0.0742138287818073 | 0.115649991879242 | 0.168563999805984 | 0.233555951882637 | 0.310487513335769 | 0.398078302603856 | 0.492768744130255 | 0.588625507090279 | 0.676929905267031 | 0.753772199676962 | 0.823106046514196 | 0.892776064568763 | 0.964026386513837 |
| Chinese_loess Means | -0.00992922030088753 | -0.00685627234272808 | -0.00263644972242058 | 0.00456650135793812 | 0.0171800768542588 | 0.0379487356429578 | 0.0695048935574124 | 0.114680636092733 | 0.176234630227922 | 0.255681643811977 | 0.353044690591995 | 0.465903186889885 | 0.587480671879438 | 0.70545496017147 | 0.811956397518547 | 0.90571144015868 | 0.990020919073498 | 1.06943656618837 |
| Japanese Means | 0.474649163724578 | 0.443432703865345 | 0.410041741204845 | 0.372273786700216 | 0.329121724356353 | 0.283030955690219 | 0.240090867441392 | 0.209371090119366 | 0.201589325186606 | 0.226840139404003 | 0.291350184139259 | 0.39377817553382 | 0.519726499901562 | 0.642890482156585 | 0.747382040461932 | 0.829407324688329 | 0.892195580878682 | 0.945280306169006 |
| Lower-MiddleY Means | -0.00682865849885968 | -0.00060483620796883 | 0.00671174779516591 | 0.0168353119264841 | 0.0318371795131959 | 0.0536523053498361 | 0.0841840911894777 | 0.125268923860985 | 0.179000307645022 | 0.247353016268821 | 0.331229454779227 | 0.42919493710873 | 0.535517759646172 | 0.639543869787761 | 0.736041145928061 | 0.82613958614389 | 0.915048715480776 | 1.00384399368775 |
| South trop ch Means | -0.00259675023606447 | -0.00369446566682659 | -0.00363747351558766 | -0.000390054740314468 | 0.00904019752001561 | 0.0284916546492706 | 0.0620623767539311 | 0.113139640613902 | 0.183967610626003 | 0.275388397726442 | 0.385135385537217 | 0.507968233105058 | 0.634326208132707 | 0.750251125909705 | 0.848539337462573 | 0.930537024584384 | 1.00231564635787 | 1.06963296926003 |
| NA Means | 0.0366402042172722 | 0.0345449445801444 | 0.0333093635431334 | 0.0345629129341452 | 0.0409323268055525 | 0.0561025771005159 | 0.0840836906480487 | 0.12805515000113 | 0.189442050381472 | 0.269535827938253 | 0.368930130834347 | 0.484509427569726 | 0.606456684388843 | 0.720252238756432 | 0.817753922638379 | 0.899780585526795 | 0.971748613294788 | 1.03906371492767 |
| Southwes amaz Means | 0.207835723880993 | 0.197019933600991 | 0.185689737006914 | 0.173364610555289 | 0.159843367441753 | 0.145872982114604 | 0.134154420414965 | 0.129609533908725 | 0.138844389475622 | 0.168981654647338 | 0.225948263430072 | 0.312180448593161 | 0.422378015176177 | 0.542599962121189 | 0.664662002239541 | 0.784132607465062 | 0.897503322879635 | 1.00679099727573 |
| C/S_Andes Means | 0.0291205853249027 | 0.0314816768181023 | 0.0344713731256695 | 0.0390464429788137 | 0.0464560633526933 | 0.0581169360685076 | 0.0760540915093959 | 0.103031413769306 | 0.142906851931112 | 0.200040121010661 | 0.27794292932455 | 0.377058478195417 | 0.491456226460239 | 0.608384004959301 | 0.721451798462705 | 0.828972489197 | 0.930712639008333 | 1.02909296975965 |
#global.SD
SD_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(SD_matrix) <- c("origin", rev(seq(4, 21, by=1)))
SD_matrix[,1] <- names(global.SD)
for(i in 1:20){
SD_matrix[i,2:19] <- global.SD[[i]]
}
kable(SD_matrix, caption= "SD values")
| origin | 21 | 20 | 19 | 18 | 17 | 16 | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| W_African_Sav SD | 0.186556185173684 | 0.159071135711959 | 0.130974576751884 | 0.101946183025733 | 0.0724951714757066 | 0.0464817364491077 | 0.0281882850407625 | 0.0212718909296576 | 0.0230484233864053 | 0.0260000666669227 | 0.0272778672869396 | 0.0272414655989988 | 0.0286145237100166 | 0.0359445310006679 | 0.0509667284519141 | 0.0718957662099432 | 0.0961865477719618 | 0.12156398251189 |
| Sudanic_Savan SD | 0.0182791317887036 | 0.0120186235797875 | 0.00857406807198337 | 0.0119997090375925 | 0.0200097740504662 | 0.0296208384871706 | 0.0390634218166876 | 0.0465027632519908 | 0.0502676895837266 | 0.0490728899403298 | 0.0422563308391379 | 0.0301155921728213 | 0.016251502147543 | 0.0230519738470708 | 0.0490762671854076 | 0.0798060762063573 | 0.111190090221826 | 0.142489654117473 |
| West Africa T SD | 0.0128551006323539 | 0.0100809743903953 | 0.0110855392306364 | 0.0157356177847031 | 0.0225339267643659 | 0.0302798766216314 | 0.0378753478987503 | 0.0441002370122807 | 0.0478360621585905 | 0.0482472418609378 | 0.045098717097735 | 0.0390673319501593 | 0.0318512113726889 | 0.0254461299782974 | 0.0263553160311764 | 0.0393349506349416 | 0.0572833097318681 | 0.0763495674818114 |
| Ethipian plat SD | 0.198659638310769 | 0.170440069868076 | 0.141737372134567 | 0.112259279537512 | 0.0864755397538964 | 0.0791421681652602 | 0.0986717402217165 | 0.129156174657549 | 0.154623472712395 | 0.166247254796878 | 0.159947337518384 | 0.135808529359533 | 0.0994473010536208 | 0.0618393023227622 | 0.0356919120363483 | 0.0377802530877869 | 0.0546872860266305 | 0.0726120006791913 |
| NA SD | 0.352588159054744 | 0.332443095173091 | 0.310979361000686 | 0.28580375866142 | 0.254547237522766 | 0.216346143796596 | 0.171980702025272 | 0.12388521679558 | 0.0763676236214709 | 0.0403854093134876 | 0.0435516738756018 | 0.0673123480707615 | 0.0842474074682211 | 0.089847619869054 | 0.0865158533690929 | 0.0824620916709309 | 0.0858640772961775 | 0.09604764429718 |
| E_North_Ameri SD | 0.00744774460623942 | 0.00688305183133896 | 0.00730077417171542 | 0.00859083907761313 | 0.0105320687647222 | 0.012800810910875 | 0.0151914913730156 | 0.0175920880323135 | 0.0199689872021909 | 0.0222421032252577 | 0.0241465903979098 | 0.0252306084926904 | 0.0249032846153157 | 0.0225138982726739 | 0.0181697023167219 | 0.0136609854949669 | 0.0121135163669742 | 0.014408338225588 |
| New_Guinea SD | 0.0890975299550152 | 0.0842301251797607 | 0.0792178489696074 | 0.0736641677047341 | 0.0671594485979282 | 0.0596008898697251 | 0.0511674521224338 | 0.0422240170397645 | 0.0331812675142621 | 0.0254299321136062 | 0.0225457007202054 | 0.0268237812941903 | 0.034037094432116 | 0.0390786810653157 | 0.0406181672896982 | 0.0401664674846312 | 0.0399927675349423 | 0.0409680658412061 |
| Mesoamerica SD | 0.0257631275045289 | 0.0187178864900728 | 0.0164812107253495 | 0.0223842305919749 | 0.0343001607587918 | 0.0492203013799705 | 0.0648668376274174 | 0.0787357939871162 | 0.0882653199553832 | 0.0913017452354271 | 0.0865857410678889 | 0.0743419528904328 | 0.0572427399139474 | 0.0406792417457969 | 0.030609493113632 | 0.0315030760258884 | 0.0386540722319049 | 0.0479903934364773 |
| N_Lowland_SA SD | 0.111921322285171 | 0.107282514497114 | 0.102185138216149 | 0.0958207981245285 | 0.0873593968043197 | 0.0763936956444485 | 0.0630414035072374 | 0.0480814037523411 | 0.0329504357663719 | 0.0219371014627934 | 0.0252549283205093 | 0.0397679138235974 | 0.0556648469342838 | 0.0682249344844884 | 0.0766047903929501 | 0.0827208213023969 | 0.0895838817981019 | 0.0974704711325054 |
| NW_Lowland_SA SD | 0.150541353866891 | 0.139669442413694 | 0.127930041006379 | 0.114205838438166 | 0.097868424572761 | 0.0797876044778232 | 0.0631128649695468 | 0.0525350681168884 | 0.0499851737421644 | 0.0520614502101073 | 0.0541664339731254 | 0.054016322802996 | 0.0513832224992821 | 0.0467970105936129 | 0.0411427020149635 | 0.0364198792817815 | 0.035220746530383 | 0.0381254553823437 |
| Sava_W_India SD | 0.00620248739705606 | 0.00409602037065882 | 0.00439841375218349 | 0.00750428764031471 | 0.012286775530684 | 0.0183044382262446 | 0.0252698645356455 | 0.0327879406646454 | 0.0403222428411792 | 0.0470625286607487 | 0.0518738127485427 | 0.0537645121764797 | 0.052235535975566 | 0.0476849083553341 | 0.0415125953320544 | 0.0349790611537332 | 0.029039548251698 | 0.0251114044895353 |
| S_India SD | 0.0085115048531172 | 0.00892906518867892 | 0.00967136329810452 | 0.0109786376784981 | 0.0132686066270033 | 0.0165765384328679 | 0.0204467930723038 | 0.0238711849145664 | 0.0256672972233344 | 0.0255092175722659 | 0.0237276529091089 | 0.0208724202116767 | 0.0181457061716282 | 0.0170301004732513 | 0.0168873351736634 | 0.0164690070990293 | 0.0154340912923534 | 0.0142507144708916 |
| Ganges_E_Indi SD | 0.00996215193394363 | 0.00691767159299303 | 0.00791433587724963 | 0.0127874855691734 | 0.0199162980428461 | 0.0287137566212716 | 0.0388527650554929 | 0.0499222314748722 | 0.0614506116687209 | 0.0728233254374634 | 0.0831495621605526 | 0.0908947825786612 | 0.0944845150273855 | 0.0928545508935936 | 0.0864932764809146 | 0.0768625361267711 | 0.065847656098237 | 0.0558662960341165 |
| Chinese_loess SD | 0.00723695180791585 | 0.00653465042035558 | 0.0067229060136461 | 0.00783510017868662 | 0.00968749896229179 | 0.0121031590893252 | 0.015015434915449 | 0.0186404465048965 | 0.0229918308268338 | 0.0276609710250012 | 0.0323906245272664 | 0.0365700364119769 | 0.0388300984431191 | 0.037141136645161 | 0.0311556397194864 | 0.0229315384793928 | 0.0156852378394716 | 0.0147092261255239 |
| Japanese SD | 0.184696194095272 | 0.167622323188328 | 0.150274927059077 | 0.131797262192788 | 0.11137082376895 | 0.0892794740360921 | 0.0676920587877228 | 0.050106916646986 | 0.0386317289017512 | 0.0331345985067326 | 0.0322763726395578 | 0.0340554941652883 | 0.0358963772505133 | 0.0361379405079019 | 0.035343175991101 | 0.0349568839559365 | 0.036562520951131 | 0.0402092257222996 |
| Lower-MiddleY SD | 0.0109221052239635 | 0.00879404210265073 | 0.0082644102269953 | 0.00985094611053262 | 0.0133492834979398 | 0.0182250521629163 | 0.0238393891749069 | 0.0293168212659387 | 0.0336239347179238 | 0.0358993675826794 | 0.0357122875762433 | 0.033336309677034 | 0.0300883638849522 | 0.0274478338873639 | 0.0261092333651014 | 0.0261920192405151 | 0.0270668753710159 | 0.0286050078015116 |
| South trop ch SD | 0.0164074485323459 | 0.013017339693911 | 0.0101269327626496 | 0.00883966436288355 | 0.0113546594405518 | 0.0176935756740366 | 0.0263781784905157 | 0.0362001386983728 | 0.0458282304601001 | 0.0540734242720755 | 0.0596617127431108 | 0.0616233193553812 | 0.0592826758862225 | 0.0525468261222113 | 0.0423612740395727 | 0.0308335216922742 | 0.0220855802940787 | 0.0229498120978104 |
| NA SD | 0.0809467917305336 | 0.0685832123328493 | 0.0557119272621779 | 0.0416794929850676 | 0.0273556018582627 | 0.0213125874557977 | 0.0336477364276112 | 0.0522922853903463 | 0.0688255532152485 | 0.0796532138108896 | 0.0829295203707789 | 0.0782121775868238 | 0.0673634830545416 | 0.0546060927712949 | 0.0437646882775398 | 0.0377464374007277 | 0.0370323136420871 | 0.0390839743471241 |
| Southwes amaz SD | 0.125774207173705 | 0.115361452592516 | 0.1042880086271 | 0.0916253149336918 | 0.0769919993552501 | 0.0614840306437265 | 0.0481765788042826 | 0.0417025596083742 | 0.0432490383179266 | 0.0477994697641187 | 0.0500604213341766 | 0.0477221942386258 | 0.0415134022061204 | 0.0341998323542236 | 0.0287859842127018 | 0.0282643985012614 | 0.0333826815416954 | 0.0417094659365358 |
| C/S_Andes SD | 0.0727773302277296 | 0.06676048268613 | 0.0609678447564098 | 0.055024169996668 | 0.0491091885343267 | 0.0446462966444592 | 0.0438008333720804 | 0.0470055485124615 | 0.0524694588100196 | 0.0582765439649155 | 0.0633094382982524 | 0.0669266582058461 | 0.068773937206034 | 0.0681584149345244 | 0.0650020480377401 | 0.0612611199316851 | 0.059915336729057 | 0.0611532400816508 |
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(j in 1:20){
plot(means_matrix[1,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[1,2:19]), type="n", xlab="year", ylab="Density", xaxt="n", ylim=c(0,1))
originI <- Pop[per.origin[[j]][, 1], ]
for(h in 1:18){
points(23- c(jitter(rep(4, length(originI[,h])), 4.5) + h), originI[,h], pch=19, cex=.1)
}
lines(means_matrix[j,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[j,2:19]), type="b", xlab="year", ylab="Density", xaxt="n")
}
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(1)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using one degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(2)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using two degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(3)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using three degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(4)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using four degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(5)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 5 degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(10)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 10 degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(12)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 12 degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(17)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 17 degree of freedom
# Load patricks productivity PCA data
load('Productivity_ALL.RDATA')
# Load origin shapefiles
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
# Extract the data
prod.origin <- extract(Productivity.ALL, origins)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, sd, na.rm = TRUE)
names(means) <- origins@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
# Plot
#pdf("productivity.pdf", 20, 30)
par(mfrow = c(5, 4), mar = c(2, 2, 2, 0))
for (i in 1:length(means)) {
plot(y = means[[i]], x = time, xlim = c(21, 4), ylim = c(ymin, ymax),
main = names(means)[i], cex.main = 1, cex.lab = 1, cex.axis = 1,
ylab = "Productivity (PCA axis)", xlab = "Thousand of years ago (k)",
pch = 20, lwd = 1, type = "l",
col = c("purple", "green")[origin.time.region[i]])
up <- sds[[i]] + means[[i]]
down <- means[[i]] - sds[[i]]
lines(up ~ time, lty = 2)
lines(down ~ time, lty = 2)
}
#dev.off()
a <- layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
layout.show(a)
frameplot <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(0, 2.25), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot_bottom <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(-0.25, 2), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot()
frameplot_bottom()
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
null device
1
###################
type_number <- 8
complex_figure <- function(type_number, i, means, sds){
if(i < 6) polygon(x=c(12,12,8.2,8.2), y=c(-1,3,3,-1), col=adjustcolor("cornflowerblue", alpha=0.4), border=NA)
if(i > 5) polygon(x=c(8.2,8.2,4.2,4.2), y=c(-1,3,3,-1), col=adjustcolor("limegreen", alpha=0.4), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- -c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq, type="l", ylim=c(-1,1))
#polygon(c(0, x_seq), c(0, y_seq), border="black", col=adjustcolor("cornflowerblue", alpha=0.5))
#abline(v= maxer - quantile(j)[2])
break_one_1 <- maxer
break_two_1 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_1, break_two_1, break_one_1, break_one_1), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.5), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 10]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 2+c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq)
#polygon(c(0, x_seq), c(2, y_seq), border="black", col=adjustcolor("limegreen", alpha=0.5))
break_one_2 <- maxer
break_two_2 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_2, break_two_2, break_one_2, break_one_2), y=c(1, 2, 2, 1), col=adjustcolor("limegreen", alpha=0.5), border=NA)
#abline(v=11)
type <- 1
if(type == 1){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
#polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
polygon(x=c(21:4,4:21), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 2){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x + 1 #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 3){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
means_long_y <- c(1,1,1,1,1, meanss)
means_long_x <- c(0:4, 4:21)
break_one <- break_one_2
break_two <- break_two_2
# polygon(x=c(break_one, break_one, 22, 22), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(1.9, 3.1, 3.1, 1.9), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
break_one <- break_one_1
break_two <- break_two_1
# polygon(x=c(break_one, break_one, 22, 22), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(-1.1, .1, .1, -1.1), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
#lines(x=c(break_one_2, break_one_2), y=c(1,3), col="white")
#lines(x=c(break_one_1, break_one_1), y=c(1,-1), col="white")
#lines(x=c(break_two_2, break_two_2), y=c(1,3), col="white")
#lines(x=c(break_two_1, break_two_1), y=c(1,-1), col="white")
#lines(4:21, meanss)
lines(21:4, meanss)
}
frameplot()
complex_figure(7, 1, global.means, global.SD)
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
axis(1)
axis(2)
quartz(width=8, height=8)
layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
par(mar=c(0,0,0,0))
# 1-4 label margins
blankplot <- function(){
plot(0,0, xlim=c(4,21), ylim=c(1, 1.25), bty="n", type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
blankplot()
blankplot()
blankplot()
blankplot()
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)
[1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" NA
[6] "E_North_Ameri" "New_Guinea" "Mesoamerica" "N_Lowland_SA" "NW_Lowland_SA"
[11] "Sava_W_India" "S_India" "Ganges_E_Indi" "Chinese_loess" "Japanese"
[16] "Lower-MiddleY" "South trop ch" NA "Southwes amaz" "C/S_Andes"
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA <NA> <NA> New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
18 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Ganges_E_Indi ... West Africa T
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
quartz
2
d <- readPNG("40962.png")
dim(d)
[1] 1000 2000 4
par(mar=c(0,0,0,0))
plot(0:360,0:360,type="n",xlim=c(20,360),ylim=c(65,295), yaxt="n", xaxt="n")
rasterImage(d, -28.5, -13.5, 388, 375, interpolate=TRUE, col=d)
axis(2, label=seq(-90, 90, length.out = 19), at=seq(1, 360, length.out = 19), las=1)
mtext("latitude", 2, line=4, at=180)
abline(h=seq(1, 360, length.out = 19), col=adjustcolor("grey10", alpha= 0.4), lwd=1)
abline(h=180, col=adjustcolor("white", alpha= .5), lwd=1)
load('PopD_all_December.rdata')
# Extract the data
prod.origin <- extract(PopD.ALL, origins_subset)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
library(matrixStats)
matrixStats v0.51.0 (2016-10-08) successfully loaded. See ?matrixStats for help.
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, colSds, na.rm = TRUE)
## new values from Bruno's GAM model (produced in script called Population_Trend_per_y.R)
means <- global.means
sds <- global.gams
names(means) <- origins_subset@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
#plot(origins)
#means[[1]] +
#sds[[1]]
#scale(as.numeric(means[[1]]), center=FALSE)
name_vector <- as.character(origins_subset@data$CONTINENT)
for(i in 1:12){
if(i > 6){frameplot()}else{frameplot_bottom()}
## customize polygons for each graph
if(i == 1){ #mesoamerica #values from Larson
complex_figure(3, i, means, sds)
}
#########
if(i == 2 ){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if( i == 3){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 4){ #Fertile crescent aka Southwest asia #values from Larson
complex_figure(8, i, means, sds)
}
#########
if(i == 5){ #loess plateau #values from Larson
complex_figure(2, i, means, sds)
}
#########
if(i == 6){ #new guinea #values from Larson
complex_figure(4, i, means, sds)
}
#########
if(i == 7){ #Eastern N.A. #values from Larson
complex_figure(5, i, means, sds)
}
#########
if(i == 8){ #Andes #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 9){ #W. African Sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 10){ #Sudanic sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 11){ #Ganges #values from Larson
complex_figure(7, i, means, sds)
}
#########
if(i == 12){ #loess #values from Larson
complex_figure(2, i, means, sds)
}
#lines(4:21, means[[i]])
abline(h = 1, col=adjustcolor("forestgreen", alpha=.5), lty=2)
# add axes to some locations
if(i == 1 | i == 7){axis(2, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
if(i == 6 | i == 12){axis(4, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
#if(i == 6 | i == 12){axis(4, at=seq(2,3, by=0.25), label=seq(0,1, by=0.25), las=1)
# axis(4, at=seq(-1,0, by=0.25), label=rev(seq(0,1, by=0.25)), las=1)
# }
if(i > 6){axis(1)} else{axis(3)}
# add text
if(i < 7){polygon(x=c(-30, -30, 30, 30), y=c(-0.1, -0.5, -0.5, -0.1), col="black")
mtext(name_vector[i], 1, line=-1.2, col="white", cex=0.5)}
if(i > 6){polygon(x=c(-30, -30, 30, 30), y=c(2.1, 2.5, 2.5, 2.1), col="black")
mtext(name_vector[i], 3, line=-1.2, col="white", cex=0.5)}
# add axis labels
if(i == 1 | i == 7){mtext("scaled density potential", 2, line=4, at=1)}
if(i == 3){mtext("Thousand years before present", 3, line=3.5, at =5)}
if(i == 9){mtext("Thousand years before present", 1, line=3.5, at =5)
}
}
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:5] = 0, 1, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 0.5, 0.75, ..., 5, 7
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
Empirical CDF
Call: ecdf(maxer - match)
x[1:2] = 0, 4
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 1, 1.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
saveToPDF <- function(...) {
d = dev.copy(pdf,...)
dev.off(d)
}
saveToPNG <- function(...) {
d = dev.copy(png,...)
dev.off(d)
}
## Try them out
saveToPDF("my.pdf", height=8,width=8)
quartz
2
saveToPNG("my.png", height=8, width=8, units="in", res=300)
quartz
2
dev.off()
null device
1